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ABSTRACT 



The radial orbit instabihty drives dark matter halos toward a universal struc- 
ture. This conclusion, first noted by Huss, Jain, and Steinmetz, is explored 
in detail through a series of numerical experiments involving the collapse of an 
>■ ■ isolated halo into the non-linear regime. The role played by the radial orbit in- 

stability in generating the density profile, shape, and orbit structure is carefully 
'nI" ■ analyzed and, in all cases, the instability leads to universality independent of 

Q ■ initial conditions. New insights into the underlying physics of the radial orbit 



instability are presented. 



Ph' Subject headings: cosmology: theory — dark matter — large-scale structure of 

O ' the universe 



1. INTRODUCTION 

Spherical infall models were introduced by Gunn & Gott (1972), Gott (1975), and Gunn 
(1977) to provide a simple theoretical framework for understanding the formation of dark 
matter halos. The essence of these models is the radial collapse of successive spherical shells 
onto a primordial density perturbation. If the initial density perturbation is scale-free the 
system evolves to a self-similar state and may be followed into the non-linear regime with 
arbitrary accuracy (Fillmore & Goldreich 1984; Bertschinger 1985). 

It is widely accepted that real dark halos form through hierarchical clustering. N-body 
methods provide a direct means to study this rather complex process; the initial conditions 
are calculable from well- motivated theories of the early Universe and the dynamical equations 
involve only gravity and the assumption that dark matter is collisionless. 

Simulated halos exhibit a number of features that suggest universality in structure and 
formation history. For example, the density profiles of simulated halos across many orders 
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of magnitude in mass can be characterized by a simple two-parameter fitting formula. This 
result was first discovered by Navarro, Frenk, & White (1997) who introduced what has 
come to be known as the NFW profile: 

r/r^l + rlrs) 

The concentration parameter c associated with this profile is R^i^/rg, where normally R^ij- = 
R200, the subscript indicating the mean overdensity relative to the background. 

Similarly, Bullock et al. (2001) showed that the cumulative angular momentum distri- 
bution can be fit by a two-parameter formula to be discussed below. Moreover, Taylor & 
Navarro (2001) (see also Dehnen & McLaughlin (2005)), found that p/a^, which they argue 
is a proxy for the "phase-space density" , is well-approximated by a pure power-law in radius. 

The self-similar infall models exhibit some of the features of hierarchically formed halos 
such as gradual growth in mass and size and virialization from the inside out. But do these 
models capture the essential physics of hierarchical clustering? The naive answer is no. Many 
predictions made by the self-similar models are inconsistent with results from simulations. 
In particular, density profiles for the halos found in the radial self-similar models are well- 
approximated by pure power-laws in radius with a logarithmic slope close to —2 (but see 
Henriksen & Widrow (1999) and Lu et al. (2005) for r^^ in the outer regions). By contrast, 
the profiles of simulated halos are shallower than r~^ in the inner regions and steeper than 
r~^ in the outer regions (Frenk et al. 1988; Dubinski & Carlberg 1991; Navarro, Frenk, & 
White 1997). 

The orbit structure also presents a problem for the naive self-similar models. The 
velocity distribution of dark matter in simulated halos varies smoothly from isotropic to 
radially-biased as one moves out from the center to the virial radius. By fiat, the orbits in the 
naive self-similar models are purely radial. More general self-similarity that permits complex 
angular dependence and velocity anisotropy has been studied in Henriksen (2004). There 
it is shown that in lowest order coarse-graining, the radial dependence of these solutions is 
not changed by the angular dependence of the underlying system (see equation A-4 of that 
paper). Thus the reason for the disagreement with the simulated profiles must be sought in 
the finer grained correction terms, but the physical nature of these terms remains unclear 
(Henriksen 2006). 

We recall that hierarchical clustering involves a complex series of merger and accretion 
events which are inherently aspherical processes. Even if some of the progenitors of a given 
halo develop steep cusps through approximately spherical infall, one might well imagine 
these cusps becoming shallower through merging and subsequent relaxation. A recent paper 
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by Dehnen (2005) suggests difficulties with this picture. Based on a theorem for phase 
mixing, Dehnen argues that in the merger of two sub-halos, the steepest cusp survives. If 
by chance some of the progenitors of a halo have cusps as steep as r~^ (that is, if they form 
through radial collapse from approximately scale-free initial conditions) the final halo will 
have a correspondingly steep cusp. The corrections to the self-similar infall referred to in 
the previous paragraph can mitigate this problem, but how are they effected? 

Here we explore whether the missing ingredient in this picture is the radial orbit in- 
stability (ROI) (Huss, Jain, & Steinmetz 1999). The ROI, as the name suggests, occurs in 
systems initially biased towards radial orbits. The effect of the instability is to drive the 
shape of the system toward that of a prolate-triaxial bar, to flatten the inner density cusp, 
and to isotropize the orbits in the inner region. Indeed, the final density profile of a system 
that has undergone the ROI is very similar to the NFW profile (Huss, Jain, & Steinmetz 
1999). The ROI is a collective effect and, as with violent relaxation (Lynden-Bell 1967), 
redistributes the energy and angular momentum of the particles. The net effect is to remove 
much of the dependence on initial conditions and hence drive the system toward a universal 
structure. 

The first hints of universality in gravitational collapse were found by van Albada (1982) 
and May & van Albada (1984) who considered the collapse of an initially spherical region. 
Although their simulations did not include the effects of cosmological expansion and used 
only 5000 particles, they were able to demonstate that a wide range of initial conditions 
leads to an apparently universal state provided that the initial state is far from equilibrium. 
The mechanism appears to be violent relaxation (Lynden-Bell 1967) and the system reaches 
a final state well-described by R^^^-law (de Vaucouleurs 1948). 

The potential relevance of the ROI in cosmological structure formation was first noted 
by Aguilar & Merritt (1990) who found that a sufficiently cold (i.e., radial) initial state is 
bar-unstable. Carpintero & Muzzio (1995) carried out a series of simulations designed to 
explore the ROI in a cosmological setting. The question raised was why, in the presence of 
the ROI, are there spherical galaxies and (presumably) nearly spherical halos? As in earlier 
studies, Carpintero & Muzzio (1995) followed the evolution of an initially spherical region 
but added the effects of the Hubble expansion and cosmological fiuctuations. They found 
that the ROI is a transient effect; that systems can develop bar-like structures early on, but 
that these structures become more spherical as hierarchical clustering continues. We argue 
that those early stages where the ROI is important are nevertheless crucial in establishing 
the r~^ density profile seen in the simulations. 

The direct link between the ROI instability and the universality of dark matter halos 
was made by Huss, Jain, & Steinmetz (1999) who carried out a series of simulations of 
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both isolated halos and of properly cosmological volumes. They note that the NFW profile 
appears to be a universal feature of all their simulated halos and that in the case of spherical 
collapse, the profile is generated by the ROI. 

In this paper, we investigate the gap between the self-similar infall models and the 
cosmological simulations. To this end, we follow the collapse of the single, isolated density 
peak from the quasi-linear phase through to the formation of a galaxy-sized halo. The model 
for initial conditions is that of a smooth spherical density perturbation with initial profile 
predicted by linear perturbation theory (Bardeen et al. 1986). A series of experiments is 
designed to isolate, in turn; the effects of the ROI, of small-scale perturbations, and of tidal 
fields. We explore the connection between the ROI and halo density profiles in detail and 
consider other indications of universality such as the angular momentum profile and the 
phase-space density profile. 

The outline of the paper is as follows: In Section 2, we review essential results from 
the self-similar infall models, Dehnen's "merging theorem" for cuspy halos, and the ROI. In 
Section 3 we describe the initial conditions and numerical methods used in our experiments 
and provide an overview of our series of simulations. We present a detailed analysis of the 
simulations in Section 4 and devote Section 5 to a discussion of the ROI. Our conclusions 
are presented in Section 6. 



2. PRELIMINARIES 

2.1. Self-similar infall models 

Early attempts to understand cosmological structure formation were based on the spher- 
ical infall model (Gunn & Gott 1972; Gott 1975; Gunn 1977) in which matter from the cosmic 
background accretes slowly onto a smooth, spherically symmetric density perturbation. If 
the density profile of the initial perturbation is scale-free, that is, ii 6pi oc r~^ (e G {0, 3}) the 
system evolves to a self-similar state in the sense that its phase-space distribution function 
takes the form 

/ (r, Vr, t) = g{t)J^ ir/n{t), Vr/V.{t)) (2) 

where r^,, f^,, and Q are power-law functions of time (Fillmore & Goldreich 1984; Bertschinger 
1985). To determine the scaling relations, note that at each time t, there is a particular mass 
shell that has reached its maximum or turnaround radius and is beginning its collapse towards 
the center. The turnaround radius scales with time as rta oc t^ where S = 2/3 + 2/ (3e). The 
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self-similar solutions of Fillmore & Goldreicli (1984) and Bertschinger (1985) are derived by 
setting r\ = rta and f^, = dr\a_/dt. (See Henriksen & Widrow (1999) for a discussion of the 
self-similar infall model based on the coUisionless Boltzmann equation. This formulation 
allows for more general types of self-similarity.) 

The density profile in the self-similar state is well-approximated by a power-law, p{r) oc 
r~^, where 



,.= ('/*"*' (3) 

1^ 2 otherwise . 

In a cold dark matter dominated Universe, the primordial density field is characterized by 
gaussian random perturbations with a power-spectrum P{k) and correlation function 

ar) = j^,jp{k)e-^'-'<fk. (4) 

Consider a ua peak of this field where a = ^(0)^'^. Peebles (1984) argued that the mean 
density a distance r from the peak is given by 



6p\ u^{r) 



PI ^ 



(5) 



For a scale-free power-spectrum, P(/c) oc A;", ^(r) oc r *^"+^) and therefore /i = (3n -|- 9) / (n -|- 4) 
for n > — 1 and /i = 2 for n < — 1 (Hoffman & Shaham 1985). 

Eq. 5 is valid for a vo extremum of the primordial density field. If one adds the additional 
constraint that the region of interest is a maximum the density profile is modified (Bardeen 
et al. 1986). The implications of this modification in the context of the spherical infall model 
were explored by Ryden & Gunn (1987). 

Naive self-similarity can be extended to incorporate non-radial motions. Sikivie, Tkachev, 
& Wang (1997), Nusser (2001), Le Delliou & Henriksen (2003), and Ascasibar et al. (2004) 
include angular momentum into the infall model while preserving self-similarity and spher- 
ical symmetry. These conditions require that the distribution of particle orbits is isotropic 
in the plane tangent to the radial direction and that the specific rms angular momentum at 
turnaround is proportional to V^M^r^. Not surprisingly, angular momentum tends to soften 
the inner cusp yielding a profile that shares some of the characteristics of the NFW profile. 
(See also Henriksen (2004) who argues that flattening of a central density cusp to a core is 
a natural flnal state of relaxation processes.) 
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Along rather different lines, Ryden (1993) found self-similar solutions for the collapse 
of an axisymmetric perturbation in an expanding Universe. These models can be prolate or 
oblate and appear to have a density profile somewhat shallower than r~^ in the inner regions. 



2.2. Merging theorem for cuspy halos 

Recently, Dehnen (2005) argued that in the merger of two cuspy halos, the steepest 
cusp always survives. That is, the cusp of the remnant halo has the same logarithmic slope 
as the steepest cusp of the two progenitors. The theorem is based on the observation that 
the excess-mass function 



D if) ^ _ [F (x, v) - /] d'^d'v (6) 

JF(x,v)>/ 

always decreases under arbitrary phase-space mixing. Here, F (x, v) is the coarse-grained 
distribution function. Under a certain set of assumptions, the excess mass function for a 
single cuspy halo with p oc r~'^ as r ^ Ois oc /~2(3-7)/(6-7) g^g j ^ qq -^vi^ere Thus, shallower 
cusps are more mixed than steeper ones and merging cannot produce a cusp steeper than 
the steepest cusp of the progenitors. Moreover, D{f) is dominated by the contribution from 
the steepest cusp of the progenitors in the limit / ^ oo implying that the steepest cusp 
survives. Dehnen's theorem is supported by numerical simulations that follow the merger of 
equal-mass progenitors with various central density laws (Kazantzidis, Zentner, & Kravtsov 
2005). 



2.3. Instability 

Dehnen's theorem would seem to imply that the progenitors of present-day NFW profiles 
have r~^ density cusps since these are preserved in the final product. To understand why 
there should be this primordial cusp profile r~^ we must explore what is missing in Dehnen's 
arguments, namely collective effects such as those associated with the ROI. We find that 
the instability can lead to a spontaneous change in the density profile and orbit structure of 
the collapsing system by changing both the energies and angular momenta of the particles 
in a coherent fashion. Although Dehnen's theorem does not require these quantities to be 
conserved, it does not envisage such globally coherent action. 

The existence of the ROI was first established in the context of anisotropic spherical 
equilibrium models (Fridman & Polyachenko 1984; Merritt & Aguilar 1985). Aguilar & 
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Merritt (1990) suggested that a similar instability was operating in their collapsing systems. 
Palmer & Papaloizou (1987), following Polyachenko & Shukman (1981), gave an analytic 
proof of the existence of the instability that invokes a resonance in low angular momentum 
orbits such that the orbital precession angle of a particle in one radial period is close to vr. 
A review of the ROI can be found in Merritt (1999). 

It is interesting to note that the ROI is not the only instability that may be operat- 
ing during the process of halo formation. The self-similar solutions of Fillmore & Goldre- 
ich (1984) and Bertschinger (1985) are unstable to purely radial perturbations (Henriksen 
& Widrow 1997). The instability is driven by collective interactions between neighboring 
streams in the phase-space flow and is similar to the instability in equilibrium models found 
by Henon (1973). Once again an instability leads to local variations in the energy of the 
particles which disrupt the strict phase mixing sheets creating a coarse-grained distribution 
function. 



3. NUMERICAL METHODS 

3.1. Initial conditions 

For the spherically-averaged density profile of our initial perturbation we use the mean 
density profile of a peak of the primordial density field (Bardeen et al. 1986; Ryden & Gunn 
1987). This density profile is determined from the power spectrum P{k) which in turn 
depends on cosmological parameters. We assume an Einstein-de Sitter cold dark matter 
dominated cosmology with the power spectrurm 

P{k)=PcMk)e-'"'"/' (7) 

where / is a smoothing scale (Ryden & Gunn 1987) and Pcdm is from Bardeen et al. (1986). 
(While the current preferred cosmological paradigm has changed significantly since the work 
of Bardeen et al. (1986) - e.g. ACDM rather than Einstein-de Sitter - our conclusions hold 
irrespective of the cosmology.) The initial density profile also depends on the parameter z/ 
where ua is the height of the peak above the background. In our simulations, we assume 
/ = 100 kpc and u = 4. 

The simulation particles begin on a cubic grid and their positions are perturbed radially 
so that the density matches the global density profile {S{r)). Velocities are started in the 
Hubble flow. The small-scale perturbations included in some of our simulations are generated 
from Eq. 7 using the Zel'dovich approximation. 



3.2. N-body code 

Simulations are performed using an iV-body code (Stiff 2003) that is based on an oct- 
tree design and mulitpole expansion to compute cell-cell interactions (Delinen 2000). The 
code employs the F2 softening kernel from Dehnen (2001) and fixed time steps. The opening 
angle is dynamically computed so that the maximum force error is constant. 

The softening length e and time step At are held constant across our suite of simulations. 
Our choice of softening length is based on the criteria established in Power et al. (2003), 

AR ■ 

g ^^ (XI 



where A^vir is the number of simulation particles inside the virial radius -Rvir- For our sim- 
ulations (see below) R^i^ ~ 350 kpc and A^vir ~ 5 x 10^ yielding a softening length of 2 
kpc. 

We use a time step set by considering the maximum acceleration of the particles at the 
first force calculation: 

At = 0.1./^. (9) 

Approximately 9000 time steps are required to evolve the systems from 2; = 30 to z = 0. 



3.3. Series of simulations 

We perform fives simulations of increasing complexity. Halos (a) and (b) follow the 
collapse of a smooth, spherically symmetric density perturbation. Halo (a) is evolved with 
only radial forces while halo (b) is evolved with the full three-dimensional force law. That 
is, in halo (a), the true force is replaced by its component in the radial direction. 

In halos (c), (d), and (e), we introduce small-scale perturbations. Halo (c) is evolved 
with pure radial forces while halo (d) is evolved with the true forces. For halo (e) we include 
an external force meant to mimic the effects of tidal fields from nearby galaxies. The potential 
which generates this field is taken to be 

<l>(r,^,0) ocr^sin^^sin20. (10) 
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4. Results 

4.1. Global properties of the simulated halos 

Table 1 provides global properties of the simulated halos. The table includes the mass 
M200 contained within the virial radius R2Q0, defined to be the radius which encloses an 
average density that is 200 times the background value, pb = SH'^/SnG. Also given is the 
"alternative" spin parameter defined by Bullock et al. (2001): 

A' = -^ , (11) 

V2MVR 

where J is the total angular momentum, M is the total mass, and V is the circular velocity 
at radius R. This definition of the spin parameter is equivalent to the usual form, A = 
J\E\^^'^/GM^^'^, when measured for a singular isothermal sphere at the virial radius and is 
much easier to compute for simulated halos. In equation (11) R, M, and V are taken to be 
R200, ^200) and V200 respectively. Table 1 also gives the mean velocity anisotropy parameter 
for each halo. 



^t 



2 



^ = '-2tr '^^' 



r 



where at (fXr) is the rms tangential (radial) velocity within the virial radius. 

While the virial masses and radii are very similar across the sequence of simulations, A 
and (3 vary greatly. In the cases where the particle accelerations are constrained to be radial 
(simulations (a) and (c)), j3 ~ 1. Non-radial forces lead to a small A' even with smooth 
spherically symmetric initial conditions. Essentially, non-radial forces drive the system to 
a prolate-triaxial structure and torques between the inner virialized region and the outer 
region generate a small amount of bulk angular momentum. The effect is enhanced by more 
than an order of magnitude when perturbations are included in the initial conditions. Of 
course, the addition of an external quadrupole field, with a suitably chosen amplitude, leads 
to bulk angular momentum comparable to what is found in full cosmological simulations. 

Snapshots of the halos at various stages in their evolution are shown in Figure 1. As 
expected, halo (a) maintains nearly perfect spherical symmetry whereas (b) develops into a 
strongly prolate system. This change in shape is the result of the ROl. Likewise, halo (c), 
though beginning from slightly aspherical initial conditions, develops into a relatively smooth 
and only mildly prolate system. Halo (d), on the other hand, has significant substructure 
and is more strongly prolate. Halo (e) is very similar in appearance to halo (d). 
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4.2. Density profiles 

The density profiles for our halos are shown in Figure 2. To construct the profiles, 
we use radial bins with equal particle number (typically 500). The difference between the 
density profiles for halos with pure radial forces and those which do not is obvious - the 
steep inner cusp of the former becomes much shallower when the ROI is allowed to proceed. 
Our results agree with those of Huss, Jain, & Steinmetz (1999) in that an NFW-like density 
profile can be generated from a smooth collapse without appealing to merging. 

The steep density profiles of halos (a) and (c) are well fit by a power-law, p oc r". 
Least-square fits to the data give a = —2.022 ±0.004 for halo (a) and a = —2.050 ±0.003 for 
halo (c). These profiles are consistent with results from the semi- analytic self-similar models 
which predict, for pure radial collapse, a final density profile oi p cc r^^ with /i ~ 2. 

The other halo profiles shown in Figure (2) are fit with an NFW profile (equation 1). 
For Halo (b), the best fit profile has a concentration parameter C200 = R200/ f^s = 6.0, while 
for Halo (d) we find C200 = 5.6 provides the best fit. The concentration parameter for Halo 
(e) is C200 = 4.5. These concentrations are lower, by a factor ~ 2, than those found by 
Navarro, Frenk, & White (1997) for halos of a similar mass. Also, the density profiles in 
our halos show a sharper transition from p ~ r~^ in the outer regions to p ~ r~^ in the 
inner core. This seems to be a common feature of halos that have undergone the ROI and 
is present in other simulations not shown here. 

The density profiles of the halos that undergo the ROI rapidly and spontaneously develop 
the NFW shape. In fact, once the number of particles within the virial radius is large enough 
to calculate an accurate density profile (Figure 3, top panel), the profile is already well fit by 
NFW formula. However, the concentration parameter of the halo stays relatively constant 
(Figure 3, bottom panel) throughout the lifetime of the system even though the virial radius 
increases steadily. The implication is that the system is evolving in a self-similar fashion. 
The ROI drives a transition in the inner regions from the self-similar radial infall solutions 
to a stable state with isotropic velocities and a softer cusp. The system in this region 
is reminiscent of the steady-state self-similar solutions discussed in Evans, N. W. (1994); 
Henriksen & Widrow (1997). 

In this connection we note that Henriksen and Widrow (1995) found a steady-state 
scale-free solution in spherical symmetry with velocity isotropy that has a density profile 
oc r"^/^ with 5 > 1. 
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4.3. Angular momentum distribution 

As with the density profile, the specific angular momentum distribution also has a 
universal profile. Bullock et al. (2001) performed a high-resolution cosmological simulations 
and calculated the mass M(< j^) with angular momentum less than j^ for individual halos. 
They found M[< j^) was well fit by the simple formula, 

M(<j;) = <! ™^o+^-^ ltJ<Jmax ^^3^ 

Mtot otherwise , 

where fj, and jo are constants and jmax = Jo/ (/^ — !)• Bullock et al. (2001) also found that 
the angular momentum in their halos is well aligned; only a small fraction (about 5%) of 
the halo mass was counter- rotating (i.e., jz < 0). A similar analysis by Chen & Jing (2002) 
found that in about half of their halos, greater than 10% of the halo mass had jz < 0. 

Following Bullock et al. (2001), we calculate the angular momentum distribution in our 
halos. In Figure 4 (top panel), we show M(< jz) along with the least-squares best-fit to 
equation (13). In general, the trends seen in the halos of Bullock et al. (2001) are also 
present here; M(< jz) oc jz for jz <^ jmax, with a smooth turnover near jz ~ jmax- However, 
equation (13) does not fit this turn over well, especially for the low-A halos (b) and (d), 
where it underestimates the mass M{< jz). The angular momentum is also not well aligned 
within these two halos; about half of their mass is anti-aligned and was excluded during 
the construction of their angular momentum profiles. Halo (e), which was evolved with an 
external quadrupolar force, has a spin parameter much larger than the other two halos. Only 
a small fraction of the mass (less than 5%) in this halo is counter-rotating and the fit to 
equation 13 is significantly better than in the other cases. 

To further illustrate this point, we consider the following generalization of equation 13 

M{< J.) = Mtot J^y° ,„ (14) 

(1 + [Jz/jo)^) 

where F controls the "sharpness" of the transition from the power-law. Figure 4 (bottom 
panel) shows our halo angular momentum profiles and the best fits to this new formula with 
F = 1.52/1.32/1.14 for halos (b), (d), and (e) respectively. The implication is that the fitting 
formula from Bullock et al. (2001) is most appropriate for cosmological halos that have a 
reasonable spin parameter and angular momentum generated by an external tidal field. 
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4.4. Velocity Dispersions 

In general, the orbits of dark matter particles in simulated halos are isotropic near their 
centres and more radial further out (Colin et al. 2000; Fukushige & Makino 2001). To explore 
the connection between this trend and the ROI, we calculate velocity dispersion profiles for 
our halos. The results are shown in Figure 5. The difference between the profiles for halos 
evolved with radially constrained orbits and those that have undergone the ROI is striking. 
The latter exhibit a turnover in a,, near r^ whereas the former have a velocity dispersion 
that increases monotonically as r ^ 0. We also show the anisotropy parameter (3 (equation 
12) as a function of radius. For models (b) and (d), j3 increases from ~ at small r to ~ 0.9 
at large r with the transition occurring at r ~ r^. The trend is similar to what is found in 
cosmological simulations although there, (3 never gets much above 0.5. 

It is worth noting that the ar profiles for halos (a) and (c) are not pure power-laws in r 
contrary to what one might expect given that these models evolve from scale-invariant initial 
conditions and yield power-law density profiles. Consider, for example, a spherical halo with 
pure radial orbits and density law 

P^'^ = ^n otherwise , ^^'^ 

This system provides an idealized model for halos (a) and (c). The distribution function 
for this system is given by / {E^ j) = FqQ (E) E~^^'^5 (j^) where O is the Heaviside function 
Fq is a constant (Fridman & Polyachenko 1984). This form for the distribution function 
assumes that the potential for r < tq is ip = AnGpoVQ Inr/rQ. For this distribution function, 
0"^ = %l){r)/2 and is therefore not a power law. 



4.5. Phase-Space Density 

Taylor & Navarro (2001), and more recently Dehnen & McLaughlin (2005), investigated 
the "phase-space density" of dark matter halos in cosmological simulations. More precisely, 
the computed the ratio p/a^ as a proxy for the phase-space density. They found that this 
quantity is universally a power-law, 

4ocr-", (16) 

with a ~ 1.84 when considering the total velocity dispersion and a ~ 1.92 when just the 
radial velocity dispersion is used. 

In Figure 6, we show p/c^ for halos (a)-(d). The curve for halo (d) is closest to a pure 
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power-law. The result seems rather curious since it is halos (a) and (c) that exhibit pure 
power-law density profiles. However as pointed out in the previous section, a^ is not a power 
law and neither will be p/a^. Evidently, the ROI changes both p and a^- in just such a way 
so as to yield a power law for the ratio p/o-^. Whether this result is accidental or telling us 
something deep about the ROI remains an open question. 



5. THE RADIAL ORBIT INSTABILITY 

The most striking feature of the final halos (b), (d), and (e) is the strongly prolate shape, 
a generic feature of the ROI. To further explore the development of the ROI we calculate 
the bar-strength parameter, defined as the relative amplitude of the bisymmetric {m = 2) 
Fourier component of the mass density averaged over some region of interest. Operationally, 
this quantity is calculated through the expression (Shen & Sellwood 2004): 



•^w = Ji 



^exp{2i(pj] 



;i7) 



where (f)j is the azimuthal coordinate of the jth particle in the sum, and N is the total 
number of particles considered. The evolution of the bar strength for halos (b), (d), and (e) 
is shown in Figure 7. In the top panel, A{t) is calculated for all particles within the virial 
radius R2oo{t)- For this choice of region, A{t) is relatively constant implying that the region 
undergoing the instability scales roughly with the virial radius and that the strength of the 
bar inside the virial radius is roughly constant in time. In the bottom panel A{t) is calculated 
using all particles in the simulation and we find that A{t) oc t" where a = 2.6,2.0, and 1.6 
for halos (b), (d), and (e) respectively. Perhaps more telling is the observation that the final 
bar strength within the virial radius is the same in each of these runs. The implication is 
that so long as the initial conditions allow for roughly spherical accretion, the ROI drives 
systems to a common (universal) final state. The ROI develops more rapidly in systems that 
are initially more strongly biased to radial orbits such as halo (b). But the final state is, at 
least in the context of the bar-strength, independent of initial conditions. 

Further evidence of the dramatic effect of the ROI can be seen in the energy and angular 
momentum distributions of particles. Figure 8 shows the differential energy distribution 
dM/dE for halos (a)-(d). (The result for (e) is similar to that for (b) and (d).) The ROI 
appears to introduce a cut-off in the energy distribution at large binding energies. Binney & 
Tremaine (1987) point out the dM/dE is primarily determined by the spherically- averaged 
density profile. Our results complement this rule in that they provide an example where a 
change in dM/dE goes hand-in-hand with a change in p. 



-14- 



As previously discussed, the ROI leads to a change in the inner structure of dark matter 
halos from r~^ with radial orbits to r~^ with isotropic orbits. We note that while the former 
resembles the self-similar state discussed by Fillmore & Goldreich (1984) and Bertschinger 
(1985) the latter resembles the self-similar state found in Henriksen & Widrow (1997). More- 
over, dimensional analysis suggests that the rms angular momentum of particles in the inner 
region scales with radius as j oc ^GrM{r). This scaling is precisely the relation required to 
preserve self-similarity, as noted by Sikivie, Tkachev, & Wang (1997), Nusser (2001), Le Del- 
liou & Henriksen (2003), and Ascasibar et al. (2004). Here, we argue that this rms angular 
momentum is generated by torques due to the evolving bar-like perturbation. 

A striking confirmation of the connection between the density and angular momentum 
profiles is given in Figure 9. The figure shows j'^/rM{r) as a function of r for halos (c) 
and (d) as well as rp{r) for halo (d). Since halo (c) is run without non-radial forces, the 
angular momentum distribution is given by initial conditions. Let r^ be the initial radius of 
a particle and v±^i be its initial transverse velocity. On average, vj_^i will be independent of 
Tj. Moreover, the relation between r^ and the final average radius of a particle is given, by 
the relation 

4:71 

—r^p, = M{r). (18) 

For an r~^ density law, M(r) oc r and therefore j oc riV_L^i oc r^'^ in reasonable agreement 
with the result for halo (c). 

Turning to halo (d), we see that j'^ evolves to the value predicted by dimensional analysis 
precisely in the region of the NFW inner cusp. That is, rp{r) and j'^/rM{r) are constant 
over the same range in r. Indeed, for halo (d), these two curves look remarkably similar. 

We conclude this section with a brief discussion of the underlying physics of the ROI. 
The literature on this subject was reviewed in Merritt (1999). In spherical systems, highly 
elongated orbits are nearly closed and therefore act like a precessing ellipse or rod. In a 
system dominated by such orbits, a bar-like disturbance grows since orbits precessing slowly 
across the perturbation are trapped. 

For our simulations this argument must work for radial orbits having an initial r~^ 
density profile with the orbit distribution of halos (a) or (c), i.e., models that are unstable 
to the ROI. Recall that the azimuthal angle for a particle travelling from pericenter ri to 
apocenter r2 and back changes by 
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A0 = 2j / ^^-^ . (19) 

in (2(V^(r')-^)-jVr'2)i/2 

is a function of S and j but can be reduced to a function of r by replacing £^ and j 
with their rms values in spherical shells. The results for halos (c) and (d) are calculated 
numerically and shown in Figure 10. 

Before discussing the significance of this result we check the numerical result by calcu- 
lating A0 semi-analytically in an r"^ halo. In particular, we assume the density to be given 
by equation 15. We take the average energy to be ■?/'(r)/2 which is the exact result assuming 
purely radial orbits. (One proves this result by using the distribution function discussed 
above, but it is to be expected in a virialized system). As for the angular momentum, we 
take {j{r)) = j(ro) {r/ro) where j(ro) is calculated from the simulations. With the change 
of variables r' -^ ^Jr^r/u we have 

r^ du 

A0 = 2j/ rr^^^—A- (20) 

where rj = 2M(ro)ro/Jo- The results of this integral are also shown in Figure 10 and agree 
reasonably well with those for halo (c). 

We propose that the key result is the near constancy of A(f) at a value close to it for halo 
(c). By contrast, A0 varies considerably with radius for halo (d). Recall that halo (c) is 
unstable to the ROI whereas halo (d) has already undergone the ROI. Suppose that there is 
a perturbation along a particular axis, say, a subhalo falling in on a particular radial orbit. 
Particles at an angle 2 (A0 — tt) from this axis but over a large range in radius will precess 
into the perturbation and act in concert to enhance it. 

Together the preceding considerations make explicit the criteria given in Palmer & 
Papaloizou (1987) for the ROI. The bias towards radial orbits is evident for halo (c), but in 
addition we have shown that the resonance condition that these authors give (esentially that 
A^ = 2 (A0 — tt) is small) can hold simultaneously at many radii. The mechanism may, in 
fact, be analogous to the bar instability in rotating disks (See Binney & Tremaine (1987) 
and references therein). 

The simultaneous development of the j'^/rM{r) = constant state shown in figure (9) 
for halo (d) may also be roughly understood in terms of these ideas. Given a bar-like 
perturbation that has Mp{r) = f{r)M{r) and noting that particles within A^ of this radius 
will be most affected by the perturbation during the next two radial orbits , we can expect 
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thein to acquire an rms j proportional to the torque multiplied by the dynamical time namely 



GfM / r3 



This yields j oc f{r-)\/ Mr. Thus if / is constant, as may be expected from our initial con- 
ditions, we obtain what is observed in the simulation. However this begs the question of 
universality if /(r) were some perverse function. We suspect that universality may neverthe- 
less hold because of the rapid growth of the bar perturbation and the nearly simultaneous 
development of j^. Resonating particles will be most strongly torqued towards the end of 
their precession, and hence at late dynamical time. By then the bar contains a large fraction 
of the system mass inside the virial radius (Figure 7, top panel) and the initial /(r) is likely 
to be forgotten. 

We conclude our discussion of the ROI by considering the phase-space distribution 
function, / (x, v), for various halos. The criteria for the development of the ROI from 
Palmer & Papaloizou (1987) is that the / diverges as j -^ 0. For a spherically symmetric 
system in equilibrium, / is a function of the energy and angular momentum. Though our 
halos are aspherical and evolving, an approximate 2-integral distribution function can be 
calculated by spherical averaging and assuming that the system is in equilibrium. / is then 
given by the mass distribution in terms of the energy and angular momentum {(P M / dE dj) 
divided by the radial period, T^ (see, for example, Binney & Tremaine (1987)). In Figure 
11, we show the results for / as a function of j at various values of E. Once again, the 
distinction is clear: halo (c) shows evidence of a divergence in / at small j suggesting a 
susceptibility to the ROI while halo (d) has a distribution function that is flatter at small j 
indicating that the ROI has already operated and driven the system to a more stable state. 



6. CONCLUSIONS 

In the hierarchical clustering scenario, halos grow through the accretion of smaller sys- 
tems and, occasionally, major merger events involving comparably-sized objects. Despite the 
apparent random nature of this process, halos exhibit a number of "universal" characteristics 
across many orders of magnitude in mass. 

The accretion phase of halo formation is reminiscent of spherical, self-similar collapse. As 
illustrated in Huss, Jain, & Steinmetz (1999) and in the present work, the missing ingredient 
appears to be the ROI. In effect, this collective phenomena prevents halos from forming with 
cusps much steeper than r~^. Our analysis of the energy and angular momentum distribution 
provides strong additional evidence of the importance of the ROI in shaping the structure 
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of dark halos. 

Recently, Lu et al. (2005) considered a simple, phenomenological model for halo growth 
that reproduced many features of simulated dark halos. In the model, halo growth takes place 
in two phases, a fast accretion phase in which orbits are isotropized, and a slow accretion 
phase in which orbits remain more or less radial. The first phase leads to the r^^ cusp 
while the second produces the r"^ outer halo. In effect, the ROI instability may be the 
physical mechanism for isotropization of orbits in the first phase of their scenario. We have 
argued that the merger theorem of Dehnen (2005) together with simulations by Kazantzidis, 
Zentner, & Kravtsov (2005) suggest the need for softening of the central density cusp by some 
coherent process. Thus the ROI is a candidate for the basic collective relaxation mechanism 
that takes steep self-similar infall cusps to the inner NFW behaviour. 

The bulk angular momentum distribution of Bullock et al. (2001) appears to be a natural 
corollary of this picture in that a quadrupole field applied to the bar-like structure that 
develops from the ROI yields an M (< jz) given approximately by equation 13. Moreover 
the rms angular momentum produced by the ROI alone is given by the dimensional argument 
discussed in the text and is not in fact too dissimilar from the bulk distribution produced 
tidally, although the amplitude is much smaller 

Perhaps the most curious result is that the ROI drives the combination p/a^ toward 
a power law in r even as it drives p away from a pure power law. Whether p/a^ tell us 
something fundamental about the ROI or whether the Taylor- Navarro result is, in some 
sense, accidental, remains to be seen. 

The extent to which the ROI operates in a full cosmological simulation (let alone the 
real Universe) may be regarded still as an open question. However if the final density profiles 
have in fact to be established very early, much doubt is removed. Figure 7 suggests that the 
ROI drives isolated systems to a common final state independent of perturbations, although 
at varying rates. Less clear is the relevance of coherent effects such as the ROI during major 
mergers. This and other questions will be left for future investigations. 

It is a pleasure to thank M. Duncan for useful suggestions. This work was supported, 
in part, by the Natural Sciences and Engineering Research Council of Canada. 
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Table 1. Collapsed halo parameters. 



Halo 



Description 



M200 (IQi^Mo) i?2oo (kpc) 



A' 



/5 



a Smooth, radial 

b Smooth, unconstrained 

c Perturbations, radial 

d Perturbations, unconstrained 

6 Perturbations, external field 



2.56 


354 


2.45 


351 


2.52 


352 


2.38 


346 


2.39 


347 



0.0 



1.0 



2.8 X 10~^ 0.81 

3.1 X 10^^ 0.99 

3.8 X 10-3 0.71 

2.6 X 10-2 0.65 



EDO 




^:= p^ 







■■■■" ■ r ■* -li: ^ "■■■ '^^' ^'k^\ ■■-:■' 

^ J. >.■:■*>.■■' .^t' ^-j^ ..4: 




^ 0.0 



z = 0-3 



- ■ -,1 . fc 

-EOO D SOQ EOO HOC] EDO 200 SOB a 200 Sflfl EOO 

Fig. 1. — Snapshots of the simulations at various redshifts: from top to bottom, z = 5.5, 
1.1, 0.3, and 0. The simulations are shown projected in the x — y plane and rotated so that 
their longest axis is along the x-axis. The columns correspond to the halo labels presented 
in Table 1. 
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Fig. 2. — Density profiles of each of tlie final lialos at z = (solid lines). Also shown (dashed 
lines) is the best-fit NFW density profile for halos (b), (d), and (e) and power-law fits for 
halos (a) and (c). 
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Fig. 3. — (Top panel) The density profile of halo (d) at early (2; = 1.72, dashed line) and 
late (z = 0, solid line) times. The former is offset by 1 dex. Arrows indicate the NFW 
scale- length r^. (Bottom panel) The concentration parameter of the NFW fit to density 
profiles for the halos as they evolve in time. The solid line indicates halo (b), the dashed 
line is halo (d), and the dotted line is halo (e). 
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Fig. 4. — The angular momentum distribution for three final halos at 2; = 0. The distribution 
was calculated as in Bullock et al. (2001), and shows the total mass M(< j^) that has angular 
momentum less than j^. (Top panel) Fits to the formula of Bullock et al. (2001) (dashed 
lines). The shape parameter for each halos is /i = 1.34 for halo (b), /i = 1.27 for halo (d), 
and /i = 1.33 for halo (e). (Bottom panel) Fits to our formula (equation 14. For halo (b), 
H = 1.45 and F = 1.52; for halo (d), /i = 1.33 and F = 1.32; and for halo (e), /i = 1.37 and 
F = 1.14. 
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Fig. 5. — Velocity dispersions, in radial bins, for the final halos. Shown is the radial velocity 
dispersion dr (top panels) and the velocity anisotropy parameter [3=1 — cr'^ /(jI (bottom 
panels). The halos with constrained radial accelerations are the upper velocity dispersions 
in the top panels, and have /9 ~ 1 in the bottom panels. 
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Fig. 6. — The phase-space density of the final halos, p/o"^ (sohd fines), along witfi tfieir best 
fit slopes (dashed lines). Curves are off-set by 1 dex for clarity. Fits are as follows: Halo (a) 
-a = -1.32. Halo (b) - « = -2.04. Halo (c) - a = -1.49. Halo (d) - a = -2.15. Halo (e) 
-a = -2.04. 
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Fig. 7. — Evolution of the bar strength. The top panel shows A calculated only within 
the virial radius r2oo, while the bottom panel computes A over the entire system. Squares 
and solid lines denote Halo (b), circles and dashed lines denote Halo (d), and triangles and 
dotted lines denote Halo (e). The bottom panel also displays least-square fits to the data, 
suggesting that A oc t", with a = 2.60 ±0.03 for the solid line, a = 2.02 ±0.01 for the dashed 
line, and a = 1.63 ± 0.02 for the dotted line. 
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Fig. 8. — The differential energy distribution, dM/dE for eacli of tlie lialos. Daslied lines 
indicate those halos (a and b) which do not undergo the radial orbit instability. 
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Fig. 9. — Comparison of angular momentum and density profiles. j'^/rM{r) as a function of 
r for halo (c) (top panel) and halo (d) (middle panel). p{r)r for halo (d) (bottom panel). 
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Fig. 10. — Change in azimuthal angle for one radial period as a function of r. In general, A0 
depends on the energy and angular momentum of the particle. Here, we replace the energy 
with angular momentum with their values averaged over radial bins of radius r. Halo (c) 
- dashed line; Halo (d) - dotted line; Semi-analytic calculation for halo (c) - solid line. A 
value of A0 = tt (long-dashed line) corresponds to a closed orbit where the particle makes 
two radial oscillations per orbit. 
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Fig. 11. — Distribution function as a function of j for different values of the energy. Top 
panel - halo (c); Middle panel - halo (d); Bottom panel - distribution function recovered 
from an N-body realization of an analytic model with / = f{E), i.e., no j dependence. Line 
types correspond to different energies. In order of increasing binding energy, they are solid, 
dashed, dotted and long-dashed. 



